	PRO RATEREC,dir,ftime,Nproc
; 
leer_input,dir,N,nu,eta,eps,epse,dt,nstep,nstg,nstdat,nstk,nstp
x=2.*!pi*dindgen(N)/N & y=x

xpX = N/2 - 1 & ypX = N/4 - 1
xpO = N/2 - 1 & ypO = 3*N/4 - 1

aa=dblarr(N,N/Nproc) ; rectangulo de NxNproc
a=dblarr(N,N)
t=0.d0

apX=dblarr(ftime+1)
apXb=apX & apO=apX & apOb=apX & recflux=apX & tiempo=apX & rate=apX
delta = apX & rate_pX = apX
; For each frame
for i=0,ftime do begin
time ='000'+string(i,format='(i1)')
if i gt 9 then time ='00'+string(i,format='(i2)')
; Compute potencial a
for j=0,Nproc-1 do begin
file=dir+'a_'+time+'_00'+string(j,format='(i1)')
if j gt 9 then file=dir+'a_'+time+'_0'+string(j,format='(i2)')
openr,1,file,/f77_unformatted ; abre el archivo binario file
readu,1,t,aa
close,1
a(*,j*(N/Nproc):(j+1)*(N/Nproc)-1)=aa
endfor
tiempo(i) = t
apX(i) = a(xpX,ypX) ; guarda el valor de a en el punto x 
apO(i) = a(xpO,ypO) ; guarda el valor de a en el punto O 
endfor
apX = apX - apX(0) ; delta_a(pX,t)
apO = apO - apO(0) ; delta_a(pO,t)
recflux = apO - apX

rate(1:ftime-1)=(recflux(2:ftime)-recflux(0:ftime-2))/tiempo(1) ; rate from fin. diff.

rate_pX(1:ftime-1)=-(apX(2:ftime)-apX(0:ftime-2))/tiempo(1) ; rate from fin. diff.

;log_rate_pX(1:ftime-1)=(alog(apX(2:ftime))-alog(apX(0:ftime-2)))/tiempo(1) ; rate from fin. diff.

;CHEQUEAR TODO ESTO QUE LO HICE CON LUCA
;PLOT,TIEMPO,DERIV(TIEMPO,ALOG(RATE_PX)) ; GAMMA(T)
; Print time and reconnection rate
print,'======================================================================='
print,'Time= ',tiempo
print,' Reconnected Flux = ',recflux
;print,' Rec. Rate (a(pX)-a(pO)) = ',rate
;print,' Rec. Rate (a(pX)) = ',rate_pX
;print,' apX(t)-apX(t=0) = ',apX
;print,' apO(t)-apO(t=0) = ',apO
print,'======================================================================='

maxrate=max([max(rate),max(rate_pX)])
window,0,xs=500,ys=500
plot,tiempo,recflux,xtitle='Time',tit='Reconnected Flux',xr=[0,max(tiempo)],yr=[0,max(recflux)],xsty=1,ysty=1
oplot,tiempo,recflux,psym=2

;window,1,xs=500,ys=500
;plot,tiempo(1:ftime-1),rate(1:ftime-1),xtitle='Time',tit='Reconnection Rate',xr=[0,max(tiempo(1:ftime-1))],yr=[0,maxrate],xsty=1,ysty=1,psym=1
;oplot,tiempo(0:ftime-1),rate_pX(1:ftime-1),psym=2

;window,2,xs=500,ys=500
;plot,tiempo,deriv(tiempo,alog(abs(apX))),xtitle='Time',tit='Reconnection Rate (log)',xsty=1,ysty=1,psym=1
;plot,tiempo,alog(abs(apX)),xtitle='Time',tit='Reconnection Rate (log)',xsty=1,ysty=1,psym=1
;plot,tiempo,abs(apX),xtitle='Time',tit='Reconnection Rate (log)',yrange=[0.001,1],/ylog,xsty=1,ysty=1,psym=1
;oplot,tiempo,recflux,psym=4
;oplot,tiempo(0:ftime-1),(rate_pX(1:ftime-1)),psym=2
end

